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ABSTRACT 


Oa ' * t equations with a set of conservative and a set of non-conserva- 
tive perturbing potentials are considered. Scheifele's DS formulation of 
these equations has dependent variables similar to Delaunay*s orbital 
elements with the true anomaly as the independent variable. Efficiency 
curves of computing cost v.s. accuracy are constructed for *dams integra- 
tors of orders 2 through 15 with several correcting algorithms and for a 
Runga-Kutta integrator. Considering stability regions, choices are made 
for the optimally efficient integration modes for the DS elements. Inte- 
grating in these modes reduces computing costs for a specified accuracy. 
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1.0 IWTROPUCTION^ CONCSAISIONS 


1«1 Introduction, Summary 

This report is one study which is part of a larger NASA effort 
by many researchers to increase the accuracy and decrease the computing 
costs of generating orbits. The results will improve efficiency of 
mission desi^, operations, scientific applications and other processes 
where orbit generation is an integral psirt. 

By a "correcting algorithm" is meant a particular predict-correct 
process such as PECE which is a prediction followed by a function 
(differential equation) evaluation, a correction, auid another eveduation. 

An "integration mode" consists of an algorithm (or variable algorithm 
method) and an order (or variable order) method. 

An " orbit generatcr" consists of the formulation of the equations 
of motion (Cowell, VOP, DS, ...), an integrator, and an integration 
mode. 

An "efficiency curve" is a plot of error v.s. cost for one problem 
and one orbit generator. 

Since the cost of a particular stepsize depends on the algorithm, 
the cost of function evaluations on the fozmilation, and the cost of CPU 
time on the ccmiputer J>'stallation, there is no universally acceptable 
cost scale. All three scales are given in our figures. 

Imagine placing a set of efficiency curves, each with a different 
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algorithm (e.g. Figures 9-12 for J^). in order of increasing influence 
of the corrector to form a three dimensional graph of h v.s. error v.s. 
algorithm. Each figure is a cross section of constant "algorithm". 


ERROR, m 
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Definition: The "efficiency envelope" is the surface constructed 

by choosing the low points in the above schematic by choosing 
that order giving the smallest error for each algorithm and step- 
size. 

Choosing the low points in each cress section, we get a contour curve 
of the efficiency envelope. This process optimizes over the vsuriable 
"order" leaving the other pareuneters of algorithm, integrator, and formulation. 
Superimposing the contour curves, as in Figure 16, we project the efficiency 
envelope onto the h-error plane. Figure 17 does the same for the Jg «ind 
lunar perturbations. 

We interpret optimizing over the parameter 'klgorithnf as specifying 

«2 

an h-algorithm plane perpendicxtlar to the error axis (at 10 m in the 
schematic) and searching the lower-right (large h) side of the inter- 
section of this plane with the efficiency envelope. We thus choose the 
jilgorithm and order which allows the largest h and the smallest computing 
time, subject to the error criteria. The variable correcting and partial 
correcting algorithms (such as used in GSFC program GTDS) essentially 
refine the grid in the algorithm dimension thus allowing a choice of a 
more optimal integration mode. 

This report accomplishes this optimization for the DS formulation [I 5 ] 
and Adams fixed step integrators. A range of a factor of 10' in h gives 
a complete idea of stability and -accuracy that is possible. In one 
study j_ll} the Cowell formulation is optimized over several fixed step 
integrators. Another study ^ij shows the DS formulation compares favorably 
to others. Two relatively new integrators that should be further studied 
due to their increased accuracy and stability are the back-correcting ^2j 
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and the cyclic [5» 2l] in particular the optimized cyclic [lU]. 

The overall effort is to search for the most efficient orbit generator 
by considering the parameters 'formulation*' and 'Integrator" thus adding two 
dimensions to the geometry for a total of Analytical studies such 
as [L2 n and |l7 ] will help fill in the now four dimensional efficiency 
envelope surface by considering relationships between formulation and 
Integrator. So will numerical studies that consider different formulations 
such as [l], [4], [9], [13], [l9], [20 ]. 

Statement of the Overall Effort: 

Minimize cost as a function of formulation, integrator, integration 

mode, error 

Subject to the Constraint: error £ k. 

The approach presently taken for this enormous, nonlinear, constrained 
optimization problem is essentially a numerical mapping routine, i.e., 
choose an orbit, formulation, integrator, integration mode, stepsize, 
and then integrate to get a point on the surface. This approach is quite 
costly due to the many parameter combinations that must be studied. Using 
emperical experience and analytical studies will reduce the number of 
parameter combinations it is necessary to consider. We also recommend 
studying the feasibility of using a more efficient numericeU. optimization 
procedure [l4] to find the overall most efficient orbit generator for 
each problem. Once done the savings in applying these generators will 
be considerable. 
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1.2 Conclusions 


1. In the DS formulation (true auioma'-y) and (time element) are 
the only variables termed "fast" and is 10 faster than 

With Jg the variation in and er^ increases to the order c* 

2 

With and lunar perturbations er^» and become 10 to 

3 

10 faister than their fastest variation is near apogee, Adams 

o 7 

efficiency decreases by about 10 ( R>K by 10 ) and stability problems 
enter. For stability, atepsize control is advisable in this vicinity, 
but this shoxild be done carefully. 

2. The dynamic and numerical stability of the OB formulation depends 

on that of the DS and on the second partials of the DS-OB generating 
function. 

3. The order 7 R*K integrator has the slope of an order 10 method, while 
Adams generally act like they are lower order than they Ideally are. 
Despite this, Adams fixed step is overall more efficient than R-K 
fixed or variable step. (At higher acc\iracies with lunar perturbations 
R*K variable step is slightly better.) For best efficiency, Adams 
variable step should probably be used near perturbations other than 
those due to the central body. 

The remaining conclusions are concerned with Adauns integrators. 

U. Numerically DS is very stable with order 1^ is stable at only 7 
steps per revolution. With the moon order 6 remains stable at this 
large h but higher orders have smaller stability regions. The stability 
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region becomes smaller with more corrections, also. 

5 . With Jg* higher orders and more correcting generally yield steeper 

curves; however, orders greater than 12 do not significeuitly improve 

efficiency. For small h all ore rs are parallel with slope 1.0 or 

2.0 depending only on algorithm. With cuid moon at all h, the 

curves for different orders generally stick to a tight bundle with 

slope 1.^ to 2.3 depending only on algorithm. This behavior is not 

typical tnuication or roundoff. The explanation is hypothesized to 

lie in the truncation error coupling of formulation to integrator. 

2 

6. With Jg, PE(CE) order 12 or 15 is most efficient for errors s 1 m. 

The most efficient for larger errors is PEC order 12 or 15 . With 
2 

Jg and moon, PE(CE) order 8 is most efficient and stable for errors 
S 100 m. For larger errors PE<' order 6 is best 
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2,0 PROBLEM FORMULATION 


2,1 Problem Sets 

The greatest use of this formulation is likely to be on highly 

eccentric orbits. The orbits considered have initial conditions: 

X(l) = .9OOOOOOOOOOOOOOD 04 «(1) = -.1120590941151220 01 

X(2) = -.lOOOOOOOOOOOOOOD 04 o(2) = -.r'>*205bl3794442D 01 

X(3) = -.80OOOOOOOOOOOOOD 03 o(3) = .2765928190018920 01 

X(4) = .2600000000000000 05 »(4) = -.728599'*650403270 C5 

X(5) = -.4123110000000000 01 ff(5) = .2418154059291120 06 

X(6) = .8000000000000000 01 0'(6) * .7163419425915420 05 

X(7) = -.2000000000000000 01 «(7) = .67876890OOOOOOOO 05 

X(8) = .1358707695979750 01 o-Ce) = .1358707095979750 01 

a = 1.5 X 10^ km, e = .96, i = l8°40 , p = .58 x 10^ secs = 6.8 days, 

where X. ^ , is the Cartesian position in km, x*. ^ „ is the 

Cartesian velocity in km/sec, Xj^ is the physical time, Xg is the total 

energy, o is the OS vector described later, and a, e, i, and 

p have the usual meanings. 

We consider two sets of perturbations: (l) the full J^, conseivative, 

term of the Earth's potentials - ^ “ (r^) | where ( - 2.627622224 x 10^®. 

(2) The above term plus the non-conservative (time varying) force due 

to the fictitious moon moving in a circular orbit at a distance 3.844x10^ km 

_6 

from the Earth with an angular velocity 2.66 x 10 radians/sec 
(p s 27.4 days) and gravitational constant times mass = 4902.66. This 
fictitious moon will have the same effect on the Integrators as the real 
moon and will avoid ephemeris difficulties in this testing phase of the 
work. 


A sketch of the geometry during the first orbit is shown in Figure 1. 
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Figure 2 graphs the Garth-satellite distance, r, and the physical tine, 
v.s. the true anomaly, T, for only. Figure 3 shows the saiie for 
Jg and the moon. The final tine specified on input is t^ ■ 1.^9*^ s 10^ secs. 
Just before three revolutions and Just after apogee. With only 

> 16.899 radians, r « 2.796 x 10^ tas and with end the moon « 16.899 
radians, r « 2.736 x 10^ km. 

2.2 Formulations 

In integrating elliptical orbits much overhead and interpolation 
error can result froii using a variable step algorithm to reduce local 
truncation errors nesir pericenter. One mauis of avoiding this is to 
use automatic time step regulation ^lo]], ^16^. Transform the independent 
variable from t to s by using 
(2.2-1) dt/ds » cr” 

Taking equal steps of size As will automatically yield smaller time 
steps, At 'cr”4s, when r is small. Thus, the choice of n affects the 
efficiency of the integration [121. 

An interesting aneJLysis by Velez [171 shows the effect on the dynamic 
stability of the transformed equations of motion and, since is 

coupled to the numerical stability of the integrator, the stability of 
the whole orbit generation process is also affected by the choice of n. 

He finds n « 1.9 best for the time regularized Cowell formulation with 
a fixed step Adams integrator. The choice is difficult since small 
n yield less dynamic stability and larger n yield less numerical stability. 

Another problem with time regularization is that to get the position 
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at a final physical time the tinffi equation (2.2-1) must be integrated 
giving timing errors. Multiplying this error by the velocity we get 
an implicit error in the position which must be added to the explicit 
position error. Several researchers, including Baumgeo:i;e, Stiefel 
and Ifecozy, have is^urijved timing accuracy and stability by adding control 
terms or energy surface considerations to the equations of motion. 

In integrating orbits that are basically two body, advantage can 
be taken of the two body aneOytic solution by transforming the dependent 
variables to the two body elements. These "variation of parameter" 
methods exactly solve the two body portion of the motion ^l6j and allow 
larger integraticxi steps to be taken. Variable stepsize is still advisable 
sin<% the elements change fhstest near pericenter. 

For this reason, authors have recently coDft>ined the above independent 
and dependent variable transformations in the KS, DS auid other formulations 
fl^J. These formulations exactly solve the two body portion and have 
automatic stepsize and energy control by including time and energy elements 
in the canonical formulation, eis opposed to adding the time equation and 
control terms onto the Cowell formuJation. Uie remaining problem with 
these formulations is that some of the dependent variables are slow and 
scxne are fast. Traditionally, "slow" means the variable is constant in 
unperturbed motion. Ilie integration step must be smaller to account for 
the fast variables. Axialytical fl^} cmd numerical averaging techniques 
are being formulated which smooth out these fast variations. 

Two Scheifele ^I^JDS formulations are considered in this report. 
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In the first) we begin with the canonical equations of motion for the 
Cartesian phase space dependent variables with time as the independent 
variable. Riase space i.. extended by adding the energy-time conjugate 
dependent variables. This allows treuis format ion to an independent 
variable whi<d> is linear in the true anomaly. Canonical transformations 
are made to the spherical, phase space then to elements similar to 
Delaunay's. Of these dependent variables, seven are slow and one is 
fast (the true anomaly). The fast variable is very fast and the equations 
of motion have periodic terms multiplied by the fast variable for non- 
conservative forces (mixed secular terms). These stability and accuracy 
reducing problems sure eliminated by further transforming the depend .t 
and independent variables. 

The DS independent variable is thus the true amomaly pxas a constant 
and tiie dependent variables are: 


or^ = 0 = true anomaly 

- argument of pericenter 
= h e longitude of ascending node 
s f s time element 



Og B G <= singular momentum 

Oy B H B third component of angular momentum 

Og = L B - energy 


fsist 

slow 

slow 

fast 

slow 

slow 

slow 

slow 


The angular momentum element (integral in Kepler's second law) is 
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used both as a dependent vatriable and in trams forming the independent 

. V 2 ds 

vauriable since from (2.2-1) r ^ relations between the 

Cartesian extended pheme space vector (2.1) amd the DS elements are: 


0 ■ sii^ (position* velocity) arccos 

I h » 

\r sin if ^ 

i-°i) 


g s arcsin 
h aurctetn 
£ = t ^ 


rl - e sin 


(2L) 

» «= G + 2r^V + 




G = llposition x velocityl' 


H * G, 


L = Xq= velocity + ~ - V 


where : 

(Gi» Gg, Gj) = position x velocity 
V B perturbing potential, 
u B gravitationail const, x (M -f m) 
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sin i B sign (H) 




p B i (G - ♦ ♦ 


II 


e 




r B p(i - e cos 0 ) 


-1 


B B 2 euKitaii 




Details of transforming the differential equations to get ^ are given 
in [16^. This is integrated to get an or value which is transformed back 
to X. 

Figures 4 and ^ graph these elements v.s. true anomaly for the first 
revolution for only, and Figures 6 and 7 for J 2 and the moon. Except 
for or^ and they are periodic auid the range of variation of or^^ is by 
far the greatest. Sven though or^ and or^ sire constauit without perturbations, 
as soon as is added their vsuriation is of the order of or^. Any emsdysis 
based on their vsu'iation remaining small must be modified, and 
also vary but quite slowly, ctj remains constant with J 2 and Og with sill 
conservative perturbations, as expected. 

2 

VIhen the moon is added, the variation of o^, er^, smd Uj is 10 to 
io3 greater thsui that of but it is still less thsui Oj^. The remaining 
elements vary slowly. The greatest variation is near T » 4.4 which is 
at apogee (the point nearest the moon's orbit (Figure 3 ) not the point 
nearest the moon (Figure 1)). 
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Since the DS foroiulation was geared toward the centred body being 
the main potential source, it is not surprising that the moon ceuises 
such gross "misbehavior” which can cause stability and truncation pro- 
blems. This points to the necessity of using stepsize control when 
perturbations other than the central body are present. 

In the second formulation considered, Scheifele has tried to slow 
the fast variables by analytically averaging the high frequency terms 
due to out of the eq]uations by applying the von Zeipel algorithm. 
Ideally, this will allow larger stepsizes near the central body irtiich 
will more than cooqpensate for the increased confuting necessary to make 
the extra transformation. Using the generating function S^, a transfor- 
mation is made from or to or', the 'bblateness elements" OB. These" inter- 

2 

me<liate"elements are almost canonical since terms of order f have been 
ignored in the transformation. 

If q = c- ^« + — — and 

(2.2-1) Sj . (- i g . i ,i„ |il . I g . 2*) 

* f pi 

then the relation between O' and or* is: 

/?^s as as as \ 

(2.2-2) or'(<y) « or + , ..., ^ ..., - gjj-j = O' + cS(or) . 
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In transforming the differential equations, we first compute ^ as for 
the DS elements. 

Flrom {2.2~2), we compute the Jacobian 


(2.2-3) 

and form 




I + 


€ n 


dty' 

ds 


30f* dtt 

ds 

dor 


ds — ds 
dor 


The dynamic and numerical stability of the OB elements thus depends on 

that of the 06 elements and on the second partials of the generating 

Ay* 

function Sj^. Integrating yields a value of «*, To find or from or', 
we have or(o'’) = - rS(o') which is a nonlinear transformation and must 

be iterated, thus increasing the computing time. This transform recovers 

the Jg high frequency terms. The transform from o to X recovers 
Uie two body portion and yields the Cartesian state. These back trans- 
forms must be done at each step in order to evaluate the forces (Sectitwt 
3.2). 


3J) EQ,UATI0M IMTEGMTIOM 


3.1 Integrators 

The orbits are generated with two basic types of integrators. The 
first is the seventh order Runga-Kutta method RK7(8) of Fehlberg [6j. 
This is the highest order R-K method in use today. An estimate to the 
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loceQ. truncation error is obtained by evaluating its leading term. 
Thirteen function evaluations per step are required. Following these 
evaluations the local truncation error estimate is compared to a user 
supplied tolerance » TQL. If the estimate exceeds TOL, a smaller step- 
size is chosen eund 12 of the evaluations are repeated with this new 
stepsize. The stepsize is reduced until the estimate is less than TOL. 
Then this value for the solution is accepted at this step and we proceed 
to the next step with a specified stepsize. An option is available to 
override the stepsize control and integrate with a fixed stepsize. 

The second integrator used is the Adams, nonsummed, ordinate form 
[8^. In this study, we \ised fixed stepsizes; the smallest being 10*'^ 
of the largest which is a greater range than conqseurable studies. 

This was done to get a complete idea of the range of stability and 
accuracy that is possible with the DS elements. The orders used rsmge 
from 2 to 15. Order 15 is three or four orders higher than 
commonly used in orbit generation. The fixed correcting 
algorithms used sure PE, PEC, PECE, and PE(CE) which require one to 
three function evaluations per step. The coefficients were computed 
in 25D, using a CDC 3350. 

The progrsun will accept cyclic integrators [5], [I8], [2l] and a 
thorough study using these should be done. The optimized ones flU*| are 
especieiUy more stable euid accurate thsui the traditionsd multistep 
integrators and relatively improve with high eccentricities. With 
stepsize control they will probably require fewer stepsize changes. 
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MAIN : Initializes Caurtesian vector, transforms to initial DS vector, 

calls either RK or DiTEG, then tramsforms finad DS to Caurtesian 
coordinates . 

DSTOCO, COrODS : Transformations between DS auid Cartesiaui vectors 

(Section 2.2). 

FORCES : Conputes potentials and forces in Cartesiaui coordinates. 

ELDEEft : Converts potentiads auid forces from Cartesian coordinates to 

be used in the DS differentiad equations for 

as 

RK78: Runga-Kutta integration from one step to the next (Section 3.1). 
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IHTEG : Uses RK78 to compute starting \raJ.ues for trculltional or cyclic 
multistep, predictor-correcter integration (Section 3*1) step- 


ping along in true anomaly (T) until Xj^ > final time (t^). Now 
consider Xj^ to be a function of true anomaly X^^ (T). We wish 
to find that T value, T^, where Xj^(T^) - t^ = 0. Using Newton's 
method: 


(3.2-1) 



- tf) 


idiere the denominator is approximated by 


X^(Tj) - 


T - T 
^i i-1 


^i+1 converges to T^.. Iteration is stopped iidjen 

~ ^1 ^ UITOL which is specified on Input. Evaluating the 
OS vector at bs will yield Xj^ « t^ euid this final solution 
will automat icadly include any timing error (Section 2.2). 

BK: Controls BK78 as it steps along in true anomaly (Section 3*1) until 

then proceeds as in INTEG. 
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The second program is for the OB elements. 



Initializes Cartesieui vector, transforms to initial DS and OS 
vectors, calls either RK or IHTEG, then transforms final OB to 
DS to Cartesian coordinates. 

OBTQDS ; Transformations between DS and OB vectors (Section 2.2). 

MATRIX : Computes the Jacobian ~ (Section 2.3). 

OBDEEQ: Computes 4^ from 4St 

— to 


le 
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INTEG and RK now integrate the OB vector. At this writing, several 
bugs .ire being worked out of the OB version of the program, so we 
have no numerical results yet to report. A preliminary version of this 
program with the Runga-Kutta integrator is due to G. Schc.'fele ^15]. 

4.0 EFFICIEMCY 


The reference solutions were generated by converting the program 
to a CDC 3330 with 2^D and running it with the R-K variable step with 
TOL B DTTOL = 10 (Section 3)* Having both Cartesian and 1)S solutions 
allows an analysis of the effects of the transformations. We are con* 
fident the reference solutions are accurate to I6D. Figures 2 through 7 

graph the solution. 

Comparison runs for the OS elements were made on a Xerox Sigma 7 
which has the same word structure as t)ie IBM 360 (l60), allowing compar* 
isons with runs on that machine. A total of 48 curves were generated 
(Figures 8 through 15 and Table 1) each with approximately 6 stepsizes, 
ranging over a factor of 10 , for a total of 288 successful runs. 
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TABLE 1 


SUMMARY OF EFFICIENCY CURVES FOR THE OS FORMULATION 


Fig. 

No. 

Perturbations 

Integrator 

Algorithm 

Orders 

8 

Jg only 

RK 

fixed step 
variable step 

7 


and moon 

RK 

fixed step 
variable step 

7 

9 

Jg only 

Adams 

PE 

6, 15 


Jg and moon 

Adams 

PE 

6, 15 

10 

Jg only 

Adams 

PEC 

2-15 

11 

Jg only 

Adams 

PECE 

2 - 15 

12 

Jg only 

Adams 

PE(CE)^ 

2 - 15 

13 

Jg and noon 

Adams 

PEC 

2 - 15 

14 

Jg and moon 

Adams 

PECE 

2 - 15 

15 

Jg and moon 

Adams 

PE(CE)^ 

2 - 15 

16 

Jg only 

Adams 

all 

envelope 

17 

Jg emd moon 

Adams 

all 

envelope 


The first vertical scale is the Eiiclidean norm of the Cartesian 

position error in meters. This includes the implicit error due to the 

timing error. The second scale is the relative error which is the error 

-4 

divided by the fineJ. r value. A relative error of 10 , for example, 

means there are four significant digits of accuracy in the satellites 
position. 
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The first horizontal scale is the "effective stepsize" h, in radiauas. 
h is normalized to be the actual stepsize for the standard PECE (2 function 
evaluations per step) and is adjusted for other correcting algorithms 
according to the number of function evaluations. This enables a direct 
comparison of computing times for different e^^orithms. The fonmila is: 


2 (arc length) 

[total # fctn. evaluations] 


which is twice the average distance between evaluations. 1%e second 
auid third scales are the tctal number of evaluations and the Sigma 7 
CPU seconds not including factors for the amount of core used or swaps. 

This will clLLow compeurisons with other machines. 

Let p = the order of the integrator and d depend on the integrator 
and on the high derivatives of the differenticd. equations thus on the 
formulation used. In the log-log graphs in the truncation limited region, 
ideadly: 

(4.0-2) log (error) ^ p log h + log d . 

The curves should be straight lines with slope = order. A good ccmibination 
of formulation and integrator will yield a siualler d, thus yielding an 
overall efficiency improvement. This is one reason for optimizing inte- 
grators l8, 193 addition to improving stability, and for considering 

the combined effects of formulation and integrator [l2, ItJ. 
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U,1 Runga - Ktttta 

In Figure 8, the points for h > .2 should be Ignored since aost of 
the function evaluations sure involved in the finishing procedure (Section 
3«2) and do not accurately represent efficiency. Otherwise, the curves 
are generally linear with both fixed step slopes ^ 10 which is three 
orders higdi^r than expected. The variable step curves are even steeper, 
so Btuch so that the choice of TQL is very critical to the propogated 
error. 

7 

The looon's perturbaticm increases the error by a factor of 10 to 

Q 

10 , possibly caused by a larger value of d (U.0-2). This accuracy 
degradation is not as great with Adams integrators. The coqmting time 
is simultaneously increased due to the additional coaqmtations . Notice 
the two CPU scales. 

4.2 Adams, Ccaaservative Perturbatiwi 

Figures 9-12 show quite a regular behavior £uid all modes remain 

2 

stable. Figure 12, PE(CE) , curves exhibit idiat mi^t be temed ideal 
behavior. VAiat we will ce^J. a "pivot point" occurs at (b, error) 

(.3> 10^ m) where all curves cross almost simultaneously. At smaller 
h, bi^er orders are more efficient; at larger, they are less efficient. 
The lines should cross due to the expected different slopes (4.0-2) but 
this does not explain why they all cross at nearly the same point. 

At the largest h run, we are taking only about seven steps for a 
whole revolution. The stability is remarkable. At even larger h, the 
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hi^er orders will become unstable first. 


As h decreases, the error reaches a minimum then increases as 
roundoff takes over causing a”diiT near (h, error) = 13«3 significant 

digits) which is almost all we can expect from this 16D machine. The 
orders 2 and 4 truncation errors are still too large for the total error 
to be seriously effected by roundoff. At even smaller h, these ^ould 
8lLso show a dip. 

Figures 9, 10, 11 resemble one another. They have idiat look like 
truncation limited regions with the pivot point and dip region moving 
to the left in h as the corrector influence increases in a definite 
trend. 


algorithm 

PE 

PEC 

PECE 

PE(CE) 

pivot h 


1.6 

.6 

.3 

dip h 

.6 

.2 

.1 

< .06 


However, these dips occiur at less than 9D> implying their cause is not 
roundoff. 

As h decreases, all orders become parallel with slope - 1 for PE 
and slope = 2 for PEC and PECB (Table 2). This is not a truncation 
limited region since the slopes are equal auad not a roundoff limited 
region since errors decrease with h. We hypothesize the explanation for 
these phenomena lies in the truncation error coupling of the formulation 
to the integrator. 
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4.3 Adams, Hon -Conservative Perturbation 

Figures 9, 13, 14, 15 with the lunar perturbation aure no^ as reguleir 

as with J2 only and the higher orders do become unstable. In Figure Vy, 

2 

the PE(CE) curves sire nearly strsd^t lines at small to medium h with 
bilker order methods being steeper as expected in a truncation region 
(4.0-2); however, the nesur equsility of the slopes wsis not expected. 

At the smsOJ-est h run, truncation is still too large for roundoff to give 
a dip. 

There is a pivot near h = .05; however, to the right there sure 
convex cmd concave curves, so the picture is not as clear as with only. 

Orders 10, 12, auid 15 became unstable sind the higher orders have 
smaller stability regions. It is surprising the order 8 remained stable 
at h giving only 7 steps per revolution. 

Figures 9> 13> l4 resonble one another. Orders 12 and I 5 always 
become unstable after having a dip at error values too laxge to be caused 
by roundoff. The point of instability (auid the dip) moves to the left 
with increaising corrections (ais with J 2 only), so the stability region 
becomes smaller with increaising corrections. 

There appear* to be no tnwction nor roundoff regions at all. The 
cxirves Just twist auround each other in a tight bxuidle with amy order 4-8 
being as efficient as euiother (2 is also close), except the 12 auid I5 
8 U% more efficient at the dip and less when they are unstable. The 
bundle hats slope approximately 1.6 for PE auid 2.3 for PEC auid PECE which 
is ainalogous »o the small h result for J 2 only (Table 2). 
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TABLE 2 


SLOPES OF EFFICIENCY CURVES 


Algorithm 

Stepsize 


Orders of Adams 

Methods 


1 

j 

! 

t 

Range 

2 

4 

6 

8 

10 

12 

15 

! 

1 



1 

2 Only 





^ 

PE 

small 



1.0 

• 


• 

1.0 


f 

medium 

- 

- 

* 

- 

- 

- 

* 


; 

large 

• 

- 

3.9 

- 

- 

- 

15.0 


1 

i PEC 

ftBfiftl 1 

2.0 

2.0 

2.0 

2.0 

• 

2.0 

1.9 


r 

t 

medium 

2.0 

2.0 

2.0 

♦ 

- 


* 


• 

large 

2.7 

^.3 

5.0 

6.6 

- 

8.0 

12.0 


' PECE 

small 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 


. 

medium 

1.8 

* 


* 



* 


. 

large 

2.8 

4.0 

4.8 

6.0 

7.5 

8.0 

17.0 


1 PE(CE)^ 

small 

2.0 

4.0 

* 

- 1.0 

- 1.0 

- 1.0 

- 1.0 


t 

mediiun 

1.7 

4.0 

6.3 

* 

* 


* 


} 

1 

large 

2.8 

4.1 

5.0 

6.6 

8.6 

8.0 

12.0 


1 










1 




and Moon 





i PE 

all 

- 

- 

1.6 

- 

- 

- 

« 


i PEC 

small 

2.0 

2.6 

2.6 

2.6 


2.7 

2.5 

1 

1 

j 

medium 

1.8 

2.2 


1.7 

- 

« 

* 


f 

j 

large 

2.4 

2.3 

♦ 

2.7 

- 

es 

9 

1 

1 

1 PECE 

small 

1.9 

2.6 

2.8 

1.8 

* 

2.7 

2.5 

1 

t 

medium 

1.9 

2.5 

2.6 

2.2 

- 

6.4 



i 

large 

2.3 

2.4 

1.8 

1.9 

- 

CD 

00 


PE(CE)^ 

small 

1.9 

3.9 

6.0 

6.7 

- 

7.5 

7.5 

J 

f 

medium 

2.3 

4.2 

5.0 

;.2 

- 

6.7 

00 


1 

i 

i — 

large 

2.5 

2.1 

2.4 

2.0 

* 

e» 

00 

i 

1 


* Large curvature in this region 
- No runs 
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4.U Efficiency Envelopes, Cong>arisons 

In the R - K cxirves (Figure 8), we see that with only the fixed 
step algorithm is more efficient. This is probably because the step- 
size control is not optimal. With the moon variable step is more effic- 
ient for errors < 10 meters, but less efficient for larger acceptable 
errors. These results and Figures 2 and 3 show that stepsize control 
is not needed for only, but the stepsize at apogee is l/lO that at 
perigee with the moon indicating that stepsize control is necessary in 
this case. 

Aligning the vertical and horizontail scales of the only curves 
allows comparison of the most efficient orders for each algorithm in 
Figure 16. There is one curve for each algorithm. The PEC graph was 
shifted slightly to the left of the PE graph due to the sli^tly greater 
CPU time. Algorithms with more corrections are generally steeper and 
there is an algorithm pivot point near (h, error) = (0.l8, 1 m). These 
phenomenon resemble those of the different orders of a single algorithm; 
more corrections cuid hi^er order both lead to steeper curves (Table 2). 
The stability intervals are very large; not exceeded by our largest h. 

Using orders greater than 12 generally does not significantly 
improve efficiency. Orders 12 and 1^ dominate Figure l6. 

The most efficient integration mode for errors < 1 m and h < .18 
2 

is PE(CE) order 12 or 15. The most efficient for larger errors is PEC 
order 12 or I5. The maximum achievable accuracy is I3.5 D on the 16 D 
machine. 
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‘Rie most efficient orders for each algorithm with the lunar pertur- 
bation sketched in Figure 17. More corrections aind higher order 

lead to steeper curves and there is a pivot near (h, error) = (0.03, 100 m), 
but these phenomenon are not as clear as with only. Ifore perturbaticxis 
tend to confhse the picture. The stability intervals are still large, 
expecially for orders ^ 8. Orders 6, 12, sund 13 dominate this figure. 

Since orders 12 and 13 become unstable and in their stability regions 
they are not much of an iiig>rovefflent over lower orders, it seems safer 
to use lower orders. The most efficient and reliable integratxon mode 
for errors ^ 100 m and h ^ .03 is PE(CE) order 8. For larger errors 
it is PEC order 6. These are the same two algorithms as without the 
moon, but stability has forced us to lower orders. 

Since the stability region for Adams PEC is ideally smaller than 
PECE and even PE fl? Table l], it is surprising that with the DS 
formulation PEC is best at large stepsizes. 

Comparing Figures 8, l6, 17, we see with J2 only RK fixed step is 

2 

as efficient as Adams PE(CE) order I3 for errors < 1 m, but for larger 
allowable errors Adams PEC order 12 is most efficient. For J2 only, the 
Adams integrator is more efficient at all stepsizes. 

With J2 and lunar perturbations RK variable step is most efficient 
for errors ^ 100 m, but for larger errors Adams PEC is better. 
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Pig. i. Earth-Moon-satellite 
geometry { roughly sketched In the 
same plane j diwing first orbit; 
time marks on Jfoon and satellite 
orbits are 10®s; distance in km; 
perigee and apogee marked P and A, 
respectively 
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1m« ANCNMMLV, nK( 

2. Reference solution: Eartb-satellite distance and physical time 

ra ’’true'* anomaly and stepniasber; without moon} every second step is 


»TH» fechnicai Meocrimdun 33-TlO 
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Hrue" anomaly and Btepnxanberj with moon; every fourth step is plotted 
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Fig. U, Reference solution J2 only: for first 
orbit; quickly varying components of the DS 
elements; scale: 1 cm = 1 unit; compare to 

Fig. 2 
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Fig* 5* Reference solution J2 only: for first 
orbit; slowly varying components of the DS ele- 
ments; scale: 1 cm « 0.001 units; compare to 

Fig. 2 
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Fig. 6. Reference solution J2 and moon; for 
firrt orbit; quickly varying components of the 
DS elements; scale 1 cm = 50 units; compare to 
Fig. 3 
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Fig. 7* Reference solution J2 and moon: for 
first orbit; slowly varying components of the 
D6 elements; scale: 1 cm = O.CX)l units; 

compare to Fig. 3 
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Fig. 10. Efficiency curves J2 only Adams PEC 
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Fig. 11. Efficiency cur’.'es J2 only Adams PECE 
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Fig. 12. Efficiency curves J2 only Adams 
PE(CE)2 
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Fig. 13. Efficiency curves J2 and moon Adams 
PEC 
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Fig. lU, Efficiency curves J? and moon 
Adams PECE 
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FUNCTION EVALUATIONS 
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Fig. 15. Efficiency curves J2 and moon 
Adams PE(CE)^ 






Fig. iT- Efficiency envelopes J2 and 
moon Adams 
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